FACILITY FORM 602 


National science rou 


GP-411 


Technical Report N 
Ionosphere Radio Lak 

N 67 3935 4 _ 

(ACCESSION NUMBER) 

_ _ 

(PAGES) 

(NA'sA CR OR TMX OR AD NUMBER) 


ELECTRICAL ENGINEERING RESEj 
ENGINEERING EXPERIMEt 
UNIVERSITY OF ILL 
URBANA, ILLINOIS 




COMPUTER METHODS FOR SOLVING A DIFFUSION 


BOUNDARY VALUE PROBLEM 


by 

V. S. Gylys 
September, 1967 


Sponsored by 

National Aeronautics and Space Administration 

NsG 24 
and 

National Science Foundation 
GP-411 


Technical Report No. 32 
Ionosphere Radio Laboratory 


Electrical Engineering Research Laboratory 
Engineering Experiment Station 
University of Illinois 
Urbana, Illinois 



TABLE OF CONTENTS 


1- Introduction 1 

2. Normalization of Independent Variables 3 

3. Finite Difference Methods 5 

4. The Backward Difference Method 8 

5. A Central Difference Method Based on Three Time Levels 11 

6. The Crank-Nicolson Method 15 

7. References 19 

Appendix A Coefficients of the Differential Equation and of the 
Boundary Conditions in Terms of the Input Parameters 

Appendix B Usage of the Computing Program 


Appendix C 


Flow Charts 



1 


1. INTRODUCTION 

This report describes three finite difference methods and a computer 
program for a boundary value problem which models the density of electrons 
in the ionosphere. Let the variables H and T represent the altitude (height) 
and time, respectively; let N = N(H,T) represent the electron density in the 
ionosphere. Then this boundary value problem can be expressed in terms of 
the parabolic differential equation 

3N 3 2 N 3n 

-g* = a(H,T ) — - + b(H,T)-gg + [c(H,T) + d(H,T,N)].N + q(H,T) , (1.1) 

3h 


where 0 < H < H and T > 0, together with the following boundary conditions : 


(a) N(0,T) = 0 for T > 0 ; 


(b) 


3n 

75H 


+ kN = 0 at H = H* for T > 0 ; 


(1.1a) 

(1.1b) 


(c) N(H, 0) = g (H) for 0 < H % H* , where g^ is a known function of H 

and i = 1 or 2. (1.1c) 

For the special case to be considered in this report, the differential equa- 
tion (1.1) can be simplified to 


3N 3 2 n 3N 

= a(H)-— + b(H)-gu + e(H).N + q(H,T) 
9H 


( 1 . 2 ) 


We shall be interested in numerical solutions of this boundary value 
problem for the values of T restricted to 0 < T < T*, where T* is some time 
constant (typically, T* = 24 hours). 

Next we shall briefly outline the contents of this report. Section 2 
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discusses a normalization of the independent variables, H and T, and then 
formally restates the boundary value problem (l 0 la) - (1.2) in terms of the 
transformed variables . Section 3 begins with an explanation of the notation 
and other concepts used in the formulation of the finite difference equations 
then it briefly mentions some theoretical aspects, such as convergence, of 
the three finite difference methods to be introduced in Sections 4, 5, and 
6. Appendices A, B, C, and D describe and illustrate the IBM 7094 computer 
program. Appendix A gives the expressions needed for the computation of 


the coefficients of the differential equation and of the boundary conditions 
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2. NORMALIZATION OF INDEPENDENT VARIABLES 


For positive constants H* and T*, set 



and 


t 


T_ 

T* 


(2a) 


Since 0 < H < H* and T > 0, we have 0 < h < 1 and 0 < t (< 1 if T < T*) . The 
partial derivatives in the equation (1.1) or (1.2) can now be expressed in 
terms of the normalized independent variables, h and t, as follows: 

3N _ 8N dh 9N 1 

cSH ~ ITh dH ~~ Dh * H* , ( 2 * 2 ) 


9 2 N 9 9N 8 9N l dh 9 2 N 1 

9H 2 = = rtlii * H* ‘ dH " 3h 2 • ^2 ' 


and 


9N _ 9n dt _ 9N . i_ 
"5 t - T3T ’ dT " "5t ’ T* 


(2.3) 


(2.4) 


Consequently, the boundary value problem will now formally be written in 
the form 


9N 9 2 N , 9N 

— = A(h) 2 + B < h % + E(h)N + Q(h,t) , 
° 9h 


(2.5) 


where 0 < h < 1 and t > 0, with the boundary conditions : 


(a) N(0,t) = 0 ; 


(2.5a) 


, v 9N i 

(b) • N = 0 at h = 1 for t > 0 ; 


(2.5b) 


(c) N(h,0) = G^(h) for 0 < h < 1 , where is a given function (2.5c) 


of h and i = 1 or 2. 
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From now on when referring to the boundary value problem we shall have 
in mind the system described by the equations (2.5) through (2„5c). Expres- 
sions for computing the coefficients that appear in this system are given in 
APPENDIX A, The same appendix also contains information on the dimensional 


units used. 
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3. FINITE DIFFERENCE METHODS 

Consider the normalized independent variables h and to Let L be a posi- 
tive integer (to be provided as an input parameter) , Introduce a rectangular 
grid over the region R, 

t > 0 and 0 < h < 1 (3.1) 

in the ht-plane by subdividing R into rectangles of sides 

Ah = -i and At 
Li 

Here the time step At is a positive real number and an input quantity. Then 

the coordinates (h,t) of a representative grid point P. are 

i, n 

h = i • Ah , i = 0, 1, 2, . . . , L (3.2) 

and 

t = n • At , n = 0, 1, 2, . . . , (3.3) 

The coordinates of P. often will be written (h., t ). 

i,n i n 

In general, we shall use the subscripts (i, n) for referring to the value 
at P^ ^ of any quantity which is a function of h and t. Thus we shall write 

N for N(h , t ), the value at P. of the exact (analytic) solution of our 
i,n i J n * i,n 

boundary value problem. To cite another example, will denote the value at 

P of the coefficient A(h) of the partial differential equation; in the 
i,n 

present case, the second subscript, n, is missing because A is independent 
of t . 

Let the numbers W represent the exact solution of any one of the three 
i, n 

finite difference schemes (backward difference, central difference, or Crank- 
Nicolson method). Here we use the adjective "exact" because as is well known 
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a numerical solution obtained on a digital computer will normally contain 

round-off error and thus will differ from the numbers W. 

i,n 

The finite difference equations for the three approximation methods used 
by us are derived in the ensuing sections of this report. At this point one 
should note only the following. 

1. Each of the three finite difference methods leads to a system 

MW , = U(W ) (3.4) 

n+1 n 

of linear equations, where M is a tridiagonal matrix. In the equation (3.4), 

« 

W represents the column vector consisting of the numbers W. (i = 1, 2, . . . 

n i,n 

L) . Our computer program solves this system of equations by using a well 
known version of Guassian elimination method specially adapted to tridi- 
agonal matrices (see, for example, the reference [2], p. 104), If further 
investigation of the computational results indicated a need of different 
techniques to solve the resulting system of linear equations, our computer 
program has been so designed that any such method (an iterative technique, 
for example) could easily replace the present Gaussian elimination procedure. 
We may also mention here that the program uses double precision floating 
point arithmetic. On the IBM 7094, this gives floating point numbers with 
the mantissas 54 bits long. 

2. No theoretical treatment of the stability, compatability (also 
called consistency by some authors), and convergence properties of the finite 
difference schemes used by us has been included in this report. These 
problems are presently under investigation and will be summarized in a forth- 
coming paper. For the time being the reader may refer to any of the four 
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references listed at the end of this report: they treat these theoretical 

questions for similar (but not exactly our) boundary value problems. 
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4. THE BACKWARD DIFFERENCE METHOD 


The backward (implicit) difference analog of (2.5) is 


W. - W. W. . - 2W. . + W. _ , 

i,n+l i,n _ A l+l, n+1 i,n+l i-l,n+l ^ + 

At 1 (Ah) 2 


W w 

. r r 1+1 > n+1 i-l,n+l i E w , o 
+ V 2lAh) J + E i i, n+1 + Q i, 


n+1 


(4.1) 


Write r = At/ (£)i) . Then 


[ (At ) Q . 


Ah 


!. , + W. ] = r[ A. - B. ] W. _ 

i,n+l i,n J L i 2 l l-l, 


n+1 


+ [-2rA^ + (At) E i 


1] W. _ + r[ A. + ^ B. ] W. . 

J l, n+1 L l 2 i J l+l, n+1 


Thus for n> 0 and for i = 1,2, ...,L-1 we can write 


(4.2) 


P. W, _ + R. W. . + S W. . . = U. . , (4.3) 

i i-l,n+l l i,n+l i i+l,n+l i,n+l ’ 


where due to the boundary conditions at h = 0 we can set 


P 1 = 0 : 


(4.4a) 


the other coefficients for i < L - 1 can be written by comparing (4.2) with 
(4.3) : 


Ah 


P = r[ A. — B. ] for i =2,3,...,L-1 ; 

i i 2 l ’ ’ ’ 


(4.4b) 


R 1 = [-2r A ± + (At) E. - l] for i = 1,2, ...,L-1 ; 


(4.4c) 
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Ah 


, = r[ A. + tj- B. ] for i = 1,2, . . .,L-] 
i L 1 2 l ’ ’ ’ 


(4.4d) 


and 


«i,n + l +W i,„> for i = 1,2,. ..,1.1 . 


(4.4e) 


Obtain the coefficients of the equation 


P T W T , _ + R t W t . = U T 

L L-l,n+l L L, n+1 L, n+1 


(4.5) 


from the boundary condition 


9N 1 

■5h + 25^ N = 0 

2 

at h = 1. The central difference representation with an accuracy of 0(h ) of 
this boundary condition at t^ + ^ is 


W T . - W T , . . 

L+l, n+1 L- 1 , n+1 1 w _ 0 

2 (Ah) 2H L, n+1 

i 


or 


W - - w + W 

L+l, n+1 H L, n+1 L-l,n+l 


(4.6) 


On the other hand, the backward difference equation (4.2) for i = L is 
^L, n+1 + »L,„’ - rt \ ' f B L 1 + 


*1-^ ♦ < A ‘> \ H + r t*L + f B l ! Vl,n + 1 


Replacing W T , in (4*7) by the right-hand side of (4 fl 6) we get 

L-fl^n+l 


(4.7) 
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V +1 + W M 


] = r 


- % B i 


] 


W T . . + 

L- 1, n+1 


Ah 


+ [-2rA L + (At) E l - 1] W L>n+1 + r[A L + ^ BjM- & n+1 + W L _ 1>n+1 


]. 


Hence 


P W. 


+ R W. , = U T 


L Vl,n+1 + n L L, n+1 " “L,n+1 ' 


(4.8) 


where : 


P L = 2rA L 5 

R l = -2rA L + (At)E L 


r (Ah) r(^h) 2 

- A L " “2H^“ B L 


= - r V 2 + ^ ] + (At) [ E L - 2^ B L ] " 1 ; 


U L, n+1 “[ ^l, n+1 


W T ] 
L, n 


(4.9a) 


(4.9b) 


(4.9c) 
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5. A CENTRAL DIFFERENCE METHOD BASED ON THREE TIME LEVELS 

Since the coefficients k 9 B, and E of the differential equation (2.5) are 
independent of t and N, a central difference representation of (2*5) can be 
written in the form 

A i 2 

A. W. = — . A. (W. . + W. + W. ,) + 

t l, n 3 n i,n+l i,n i,n-l 


B. 

l 


—=■ • A (W. , + W. + W. ,) + E. . W. + Q. 

3 Ti i,n+l i,n i,n-l 1 i,n x,n 


(5.1) 


Here 


„ W. . . - 2W. . + W. _ . 

aV . = M ^ 

^ ,J (Ah ) 2 


(5 . la) 


W. . - W. . . 

aw = - x . +1 -^ izJjl , 

Vi,j 2(£h) 


(5.1b) 


and 


W. . - W. ... 

AW _ _i il+l hJzl 

A t W i,j 2 (At) 


(5.1c) 


9 2 W 

9W 

9W 

5 9 
9h- 


and — — 

9t 


at (h,t) = (h f ,t ). By expanding (5.1) in terms of these approximations and 

by, as before 5 writing r = , we get 

(Ah) 


3 (W , - W. , ) = 2A.r[ (W. . - 2W. . + W. , ) + 

i, n+1 i , n- 1 l 1 i+l,n+l i,n+l i-l,n+l 


+ (W - 2W. + W. , ) + (W. . - 2W, , + W. ■.-)] + 

i+1, n i,n i-l,n i+l,n-l i,n-l i-l,n-l 
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* B i f t (W l + l,n + l - W l-l,n + l> * < W i + l,n ' "i-l,n> + 


* ‘"i^n-l ' + "“"IVi.n + Q i,n ] 


r[2A t - OSOB,] V,,.,! + l _4rA i " 3 > W i,n*l + 


* r[2 Ai * <Ah>B.] W 1+1>Il+1 - -(2A 1 r)[tW 1+ i n + 


2(W. + W. n ) + (W + W 

i,n i,n-l i-l, n - i-l,n-J 


(B i l>' (W i + l,n + *itl,n-l> ’ '*1-1, n + *i-l,n-l ! 


[ 8 <*>«1 W i,n + - 3 W l,n-1 


Thus we obtain L-l equations 


P i W i-l,n+l * B i W i,n+l + S i W i*l,n + l = U i,n + 1 


for n > 1 and i = Since » 0 j = 0,due to the boundary condition 


at h = 0, we can set 


P ! = ° : 


the other coefficients for i < L-l are: 


P i = r[2A ± - (Ah)B.] , i = 2, 3,..., L-l 
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R. = -[4rA i + 3] , i = 1,2,..., L-l ; 


(5.4c) 


S = r[2A. + (Ah)B i ] , i = 1, 2, . . . , L-l 


(5.4d) 


and 


, = -ri[2A. + (£h)B.](W. + W. , _) + 

i,n+l V l i l+l, n l+l, n -1 

+ [2A. + (£h)B.](W. , + W. . ,) - 4A. (W. + 

L i i i-I, n i-I, n -1 i i,n 

,)l - 6 (At ) [ E . W . + Q. ] - 3W. . , 

i,n-l I 1 1 , n 1 , n 1 , n - 1 * 


+ W. 


(5 . 4e) 


i = 1, 2, „ , I>- 1 


Remark: to evaluate U, , , set W_ =Vf = 0, 

1, n+1* 0, n-1 0,n 

In order to obtain the coefficients of the L ^ equation, 


P T • W T , , + R t - W T = U _ , 

L L-l, n+1 L L, n+1 L, n+1 


(5.5) 


use the central difference representation 


W T i • “ W t 1 i 

y J } J , ™ _ o 

2 (Ah) 2 H x L, j 


of the boundary condition at h = h =1; from this difference equation we get 

L 


w. 


L+l, j 


-^W . 

L,J 


+ W. 


L-l, j 


(5.6) 


Next let i = L In "(5.2) and replace the ff by the right-hand side of (5.6) 

L+i , n 

with j = n - l,n or n + 1 as required. Thus, 
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r [ 2A. 


+ 


+ 


It follows 


and that 


L - + [ - 4M L - 31 W L,n + l + 


r[2A L ♦ f- W t;ii+1 + » L . 1>n+1 ) = 

- r|[2A L + ^(W L n A + <Vl,n * Vl,n- 1 » + 


[2A l - (Z*)B L ](W L . 1;n + W L _l,n-l 


) “ 4A (W + W 
L L, n L 


,-!>} - 


6 <“>l E L*L,n * V’ ' 3 \n-l 


that 


P L = 4rA L ’ 


Ah (At)B T 

R l =-[2rA L (2 +s -) + 3] , 


(5.7a) 


(5.7b) 


{■ 


U = - rJ 4A (W + W .) - [2A (2 + + 

L, n+1 I L Lr 1^ n L~ 1^ n 1 L 


+ B ] (W + W 

L > L, n L 


,n-l ) ] - 


■ 6 <At) I E L W I.,n + ^L.n 1 ' 3 \,n-l ' 


(5.7c) 
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6. THE CRANK- NICOLSON METHOD 

In the present case, the finite difference representation of the differ- 
ential equation (2.5) is 


W. - W. A. 

ir5 . . _i . Aw. ... ) + 

At 2 n i,n+l i,n 


B. 


E. 

l 


+ • A (W. , + W. ) + — (W. + W. ) 

2 Ti i,n+l i,n 2 i,n+l i,n 


Q. , + Q. 

i,n+l i,n 

+ 2 ' 


( 6 . 1 ) 


where the operators A^ and A^ are defined by (5.1a) and (5.1b), respectively. 

Then writing r = — — — , we have 

(Ah) 


2[W. , - W. ] = A.r[ (W. . - 2W. . + W. _ . ) + 

i,n+l i,n 1 l+l, n+1 i,n+l i-l,n+l 


B. 

+ (W. . - 2W. + W. . )] +‘5 i (^r)[(W. , . - W . .) + 

l+l, n i,n l-l, n 2 Ah i+l,n+l i-l,n+l 


(W. . - W. _ )] +[E. (W. 

l+l, n l-l, n l ’ 


. + W. ) + (Q. . + Q. 

i,n+l i,n i,n+l l.n 


) ]At 


which leads to 


2 [W. . - W. ] = A. r[ (W. , _ + W. . ) + (W. . + W. . ) 

i, n+1 i,n l l+l, n+1 i+l,n i-l,n+l i-l,n 


D • A + 

-2(W. + W. )] + -=i(^)[ (W. _ + W. ) - 

i,n+l i,n 2 Ah i+l,n+l i+l,n 
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- (W , + W. . )] +[e.(W. . + W. ) + (Q + Q )jAt . 

v i-1, n+1 l-l, n L i l , n+1 i,n i,n+l i,n 


Therefore, 


r r a ir— )]W , + [-(2 + 2rA. ) + At.E.]W. _• + 

L A i 2 W J 1-1, n+1 L i i i, n+1 




+ ' 2 - 2rA i * *V W i,n 


* + -T ( S; >IW i + l,n + <«l,n+l 


Hence for n> 0 and i = l ; 2 > ..p,Lrl we obtain L-l equations 

P lVl,n+l + R i W l,„+l * R i W l+l, n+1 = U i,n+1 ' 

Due to the boundary condition at h = 0, we have W = 0 for all j 
can set 

P ! = °- 

Then the other coefficients for i < L-l are as follows : 


p i = r t A i ■ 1 = 2 > 3 >'--> L ~ 1 ; 


( 6 . 2 ) 


(6.3) 


Thus we 


(6.4a) 


l ± = AtE - 2(1 + A^), i = 1,2,..., L-l ; 


S. = r[A^ + Og- )B^], i - 1,2, ...,L 1 ; 


(6.4b) 

(6.4c) 

(6. 4d) 


and 
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U 


i,n+l 


r [ A . ~ C^)B ]W + [AtE. + 2(1 - rA.)]W. + 

i z i l- 1 , n i i i , n 


+ r[ A. 4- <^)B. ]W. , + (Q. 4- Q. )At , 

i 2 i l+l, n n+1 1, n 9 


i — 1,2, 90 „ , L- 1 


Remark: to evaluate U , use IV = 0 e 

1, n+1* 0,n 


th 


In order to obtain the coefficients of the L equation 


P T • W T _ + R • W T , = U 

L L-l,n+l L L, n+1 L, n+1 9 


use the central difference representation 


W, , . - W. 

2 (Ah) 


J k" 1 ? J + — ^ _ 0 

ft / , ^ \ r ”, “ ^ 


2H X L,j 


(6.4e) 


(6.5) 


of the boundary condition at h = h =1; this difference equation implies ' 

L 


W , . = - <^)W. . + W. . . . 

L+1,J L,J L-1,J 


( 6 . 6 ) 


Next set i = L in (6.2) and replace the W .by the right-hand side of (6.6), 

L+l, J 

where j = n, n+1. This gives us 


[rA L " 2 t S )W L-l,n+J + f" 2(1 + rA L> + AtE L ]W L,n+l + 

B 

+ [rA + ^<t£) ]*[-<£=)», , + W_ . J = 

L 2 Ah H L, n+1 L-l,n+l J 


-[■ 


B. 


rA L - T&Kl.n + t 2<1 - rA L> * AtE L 1W L,n + 




+ Vl,J + <«L,„ + 1 * 
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It follows 


and that 


that 


P_ = 2rA 


L * 


[-2 - 2rA L 


= -[rA L (2 * §5) + + 2 - AtE L ] , 


B, 


L ,At ^ 
1 


(6.7a) 


(6.7b) 


U L, n+1 - -{ [2rA L ] Vl,n * [2 ‘ !r \* 4 'V rA L < f[ ) ' 

- + "W + ‘‘l,**'} 

" ~ {j 2 rA L^ W L-l,n * [~ rA L <2 + IT*^ + 

* 2 ♦ AtE L ]W L n * * Xn^J ■ 


(6.7c) 
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APPENDIX A: COEFFICIENTS OF THE DIFFERENTIAL EQUATION AND OF THE BOUNDARY 

CONDITIONS IN TERMS OF THE INPUT PARAMETERS 


A.l The Input Parameters and their Dimensional Units 


In order to specify the dimensions of the input parameters and the 
coefficients of the boundary value problem, let ft and ft denote the units 
of height {length) and time, respectively. The normalized dimensional units 


for h and t used by the program internally are 1\ = 1,000 kilometers and 
lft = 24 hours, respectively. The units of h and t to be used for external 
input are described in APPENDIX B. 

The following parameters must be provided as a part of input information 
H Q in ft] , 

H x in ft] , 

D 0 in ft 2 /T] , 


Q 


0 


in 


(V 


-1 

T 


> 


Nq in ^number of electrons /\ ] , 

r a_ 1 

P 0 in [t 1 , 

F = 1 or 2 

angle $ in radians , 
angle 6 in radians , 


m 


[T] 


A. 2 Coefficients of the Differential Equation 


First let us define the following quantities : 

h-H A2 ^ 

D i (h) = D 0 * EXP [“h~J in /T ^ > 
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[t ] ; 

<|> 4 cos 6 • cos t') , 

where t r is in radians and for the unit of time in r( = 24 hours), t f = k • t 
with k = 27T radians/T. 

Then, in our special case, the coefficients A, B, E, and Q of the 
differential equation 

9n 9 2 N 9n 

= A(t,h) — - + B(t,h) -gjj + E (t,h, N)N + Q(h, t) 

9h 


-F(h-H ) 

(3(h) = . EXP[ — 2_] in 

cos X = (sin 4> • sin 5) + (cos 


are defined by 


A(h,t) = A(h) = D L (h) 


• rO 2 A 

m [\ /tJ 


> 


B(h,t) = B(h) = 


3D X (h) 

~2H, 


in [ X/t ] , 


D (h) »_i 

C(h,t) = C(h) = -(3(h) + 5 in [t ] , 

2(H 1 ) 2 


D(h,t,N) = 0 , 

and so 

D (h) 

E(h,t,N) = C(h,t) + D(h,t,N) =- (3(h) + — 

2(H) 


is a function of h only. Q(h,t), the only coefficient dependent both on h 
t, is of the form 


( h-H h-H A ^ A 

1 - (sec X)EXP[ - — U in [\ 7t] when |x|< ir/2 
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= 0 when JxJ > 7 r/2 . 


A. 3 The Boundary Conditions 


The boundary condition N(h,0) = G^, where i = 1 or : 


0 < h < 1 is either 


V h) = N 0 ' EXP {l U ' sf' - 


EXP[ 


-2(h-H Q ) 

h" 


or else 


when 0 < h < 


o 2 (h) = 


N o i 


h-H 

EXP[- -gj— ] -EXP[- 




Both G i (h) represent the number of electrons/)^ o At the 
specify which G^i =1 or 2) will be used* 

For t > 0, the boundary condition at h = 1 is 


3n 1 

i + <2ir )N = 0 > 


where is an input parameter already mentioned in the 
The boundary condition at h = 0 for t > 0 is 


, at t = 0 and for 


] - 


when H Q < h < 1 


input time one should 


present appendix. 


N(0,t) 


0 
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APPENDIX B: USAGE OF THE COMPUTER PROGRAM 

B.l Computer Program 

First let us review some basic features of the computer program for the 
numerical solution of our boundary value problem which should he of interest 
to its prospective user. Almost all of them have already been mentioned in 
and are scattered through other parts of this report. However, in order to 
make APPENDIX B relatively independent from the rest of the report, we 
summarize this useful information below: 

1. The program has been coded in the symbolic language (SCATRE) for the 
IBM 7094 computer of the Digital Computer Laboratory of the University of 
Illinois. 

2. Double precision normalized floating point arithmetic is used inter- 
nally. The mantissa of a double precision floating point number on the IBM 
7094 is 54 bits long. 

3. Three finite-difference methods (backward difference, central dif- 
ference, and Crank-Nicolson) are incorporated in the program. The user must 
indicate by an input signal which one of them he intends to use. How to 
choose values for this signal is explained in the description of the input. 

4. The central difference method is of experimental nature. Its final 
acceptance or rejection will depend on the results of further testing and 
theoretical work. This method uses three time levels and thus has to be 
started by some two-times- level method which in our use happens to be the 
backward difference method. 
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B.2 Composition of the Input Data Deck 


By a run we mean an integration (numerical solution) sweep executed from 

the initial time t = 0 through some final value of t, say t . Each run may 

FN 

t h 

be segmented into one or several stages . One specifies a stage, say the i— 

by fixing a constant time step, At^, to be used throughout it and by giving 

its final (or end) time point t . In a run consisting of N stages, the t . 

Fi Fi 

must be ordered such that 


0 < t < t 
FI 


F2 


< t 


FN 


During the execution time the east stage of a run will be recognized by the 
current value of the so-called exit flag : the value of that flag associated 

with the last stage should be = 1 ; for any other stage, it should be = 0. 

A machine job will consist of one or several runs depending on for how 
many runs the input data has been furnished. 

Using the terminology introduced above we can now speak of the stage or 
run data subdecks of a job input data deck. The composition of a typical 
job input data deck is illustrated by the following diagram. 



25 


$ Data card. 


*\ 

. Data subdeek for the lat run. 


■) 

'• 


Data subdeck for the 2nd run. 


•*» *»• 


► 


Bata subdeek for the last run. 


The end-of-job card. 


Fig. B.l: Composition of a job input data deck 
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The "$ DATA” and the end-of-job cards are always of the same fixed for- 
mat and should always be present in the data deck* 

Warning: the user should become familar with the current rules of the 

Digital Computer Laboratory concerning other "$" system cards that may have 
to be added to or omitted from the job input data deck. 

The composition of a run input data subdeck is as follows: 


CARD 

Run ID (title) card: 

Run control card : 

Card 1 of floating point numbers : 
Card 2 of floating point numbers: 
Card 3 of floating point numbers: 
Stage (#1) card: 

Stage (#2) card : 

• # * • • • # « « 

Stage (#N) card : 


CONTENTS 

an arbitrary alphanumeric message in 
columns 2-72; "l" should be punched in 
col. 1. 

L; the integrati on- mode- signal ; the 
write-mode-signal; BCFLG ( = the 
boundary condition flag). 

H 0 ; H 1 ; V (V 

; 6 ; r . 

Fi N 0 ; <J 0 . 

t ; At x ; At ; SWMODE (= ‘the stage- 
write-mode signal); EXFLG (= the exit 
flag) = 0. 

t ; At 2 ; At ; SWMODE (= the stage- 
write-mode signal), EXFLG ( = the exit 
flag) = 0. 

t ; At ; At ; SWMODE (= the stage- 
write-mode signal), EXFLG ( = the exit 
flag) = 1. 


Fig. B.2: Composition of a run input data subdeck 
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The following table provides sufficiently detailed information needed to 
design and punch the input data cards. In the event that the user is not 
familiar with the definitions of the I-, E-, and H-fields as used in input 
format statements, he is advised to consult the SCATRE manual issued by the 
Digital Computer Laboratory of the University of Illinois, 


FIELD 


COLUMNS ADJUST 
& TYPE 


TYPE 

OF INPUT 
QUANTITY 


SYMBOLIC 
NAME OF 
STORAGE 
LOCATION 


INPUT QUANTITY 



Pos „ integer 


should be punched 


02-72 Arbi- 
( 71H) trary 


73-80 


lphanumeric TCAfta + 6 


Blank 


07-12 

( 16 ) 


Right 



19-24 Right Integer 

( 16 ) 


n arbitrary ID message 71 alphanu- 
eric characters long (including 
lanks); for example, 'THE ELECTRON 
ENSITY IN THE IONOSPHERE, STUDY 
OoXX". 




arameters ) 

the no, of /Sh subintervals in 
< h < 1 0 Range: 4 < L < 256 a 


The flag for determining the mode of 
the numerical solution. 

The value 0 backward difference 
method , 

The value central difference 

method * 

The value 2^ Crank-Nicolson method. 


The control flag for writing the in- 
put data and the results of the pre- 
liminary calculations „ 

The value 0=^ the input data, hj-, 

and Hj- tables written; 

The value 1 ^ in addition, other pre- 
liminary calculations 
written. 

Normally, the user should use the 
value == 0 o 


rhe control signal for the choice of 
boundary conditions at t =0. 
rhe value Os^G^Ch); 

[he value 1=^ G^(h). 


25-80 


Blank 
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FIELD 

TYPE 

COLUMNS 
& TYPE 

ADJUST 

OF INPUT 
QUANTITY 


SYMBOLIC 


C|ard 1 ofj floating pcjint numbers 


INPUT QUANTITY 


HI H 


1 in [1,000 km]. 


01-18 Right FI. point HO H in [1,000 km]. 

(E18.8) 0 

19-36 Right FI. point HI H in [1,000 km]. 

(E18.8) 1 

37-54 Right FI. point D0U D in [(1,000 km) 2 /hr] 

(E18.8) 0 

55-72 Right FI. point BETA0U (3 n in [hr" 1 ]. 

(E18.8) 0 


73-80 


Blank 


Card 2 of 

floating po 

int numbers 

01-18 

(E18.8) 

Right 

Fl. 

point 

PHI 

19-36 

(E18.8) 

Right 

FI. 

point 

DELTA 

37-54 

(E18.8) 

Right 

Fl * 

point 

TAUU 

55-80 



Blank 



Card 3 of; floating point numbers 



01-18 

(E18.8) 

Right 

Fl. 

point 

F 

19-36 

(E18.8) 

Right 

Fl. 

point 

NO 

37-54 

(E18.8) 

Right 

Fl. 

point 

QUO 

55-72 

(E18.8) 


Blank 


1 

L , ...... , 

1 



F = 1.0 or = 2.0 [dimensionless]. 

N q in [no. of electrons/(l, 000 km?]. 
Q 0 in [ (1,000 km)" 3 hr" 1 ]. 
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COLUMNS ADJUST 
& TYPE 


TYPE 

SYMBOLIC 


OF INPUT 

NAME OF 

INPUT QUANTITY 

QUANTITY 

STORAGE 

LOCATION 


Integer SWMODE The control signal for writing the 

output of the stage. 

The value - 0 or 1^ standard out- 
put ; 

The value >2 in addition, 

intermediate quantities 
dependent on t written; 
The value > 3 ^ in addition, 

1 “ the coefficients of the 

j system of linear equa- 

tions written* 

formal ly, the user should use the 
alue - 0 or L 


Integer EXFLG The flag for indicating the last 

stage of a run* It should be = 1 
for the last stage and = 0 for any 
other stage* 


02-11 

(10H) 


12-80 


TCARD + 6 phe message ? END OF J0B ? * should be 
[punched. 


Blank 
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APPENDIX C : FLOWCHARTS 


The flow charts contained in this appendix are presented in the same order 
as the corresponding subprograms appear in the program. The name by which a 
subprogram, including the master control, is referred to is its symbolic entry 
address. For the convenience of the user, we list next the names of all sub- 
programs in the order indicated above : 

MCNTRL (the master control) 

ERREX 

IAOOO 

IBOOO 

SBOOO 

EQOOO, EQ100 

MEOOO, ME100, ME200, ME300, ME400, ME500, ME600 

PROOO, PR100 

DEOOO 

WCOOO 

WTOOO, WT100, WT200 

WPOOO 

WXOOO 


Each subprogram is entered from a higher level subprogram (the master con- 
trol in this respect is the highest level subprogram) by using a basic linkage 
of the following type: 

a TSX (name of subprogram), 4 
a + 1 (normal exit) 

There are subprograms which do not return the control to the location a + 1. 
Such cases are clearly indicated, 

A reference to a lower level subprogram from a higher level subprogram is 


indicated by 
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(name of the subprogram referred to) 


brief description of the function of 
the subprogram referred to 


in the flowcharts contained in this appendix Q 



